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ABSTRACT 

Context. This paper is motivated by the recent detection of many extremely metal-deficient ([Fe/H] < —3) stars in the 
Milky Way. 

Aims. We explore numerically the chemical, thermal, and dynamical evolution of a shell formed by a high-energy su- 
pernova explosion (10^^ erg) in dwarf protogalaxies with total (dark matter plus baryonic) mass 10^ Mq at a redshift 
z — 12. We consider two initial configurations for the baryonic matter, one without rotation and the other having 
the ratio of rotational to gravitational energy /3 = 0.17. The (non-rotating) dark matter halo is described by a quasi- 
isothermal sphere. The latter choice is motivated by recently proposed mechanisms for rapid flattening of a central 
cuspy region in dwarf protogalaxies. 

Methods. We use a finite-difference numerical hydrodynamics code to simulate supernova explosions in dwarf proto- 
galaxies with axial symmetry. The advection is treated using a third-order piecewise parabolic scheme. The heating 
and cooling processes in the gas are taken into account by solving numerically the rate equations for main atomic, 
molecular and ionic species in the primordial gas. 

Results. We find that the dynamics of the shell is different in protogalaxies with and without rotation. For instance, the 
Rayleigh- Taylor instability in the shell develops faster in protogalaxies without rotation. The fraction of a blown-away 
baryonic mass is approximately twice as large in models with rotation (~ 20%) than in models without rotation. We 
argue that these differences are caused by different initial gas density profiles in non-rotating and rotating protogalax- 
ies. On the other hand, the chemical evolution of gas in protogalaxies with and without rotation is found to be similar. 
The relative number densities of molecular hydrogen and HD molecules in the cold gas (T < 10"^ K) saturate at typical 
values of 10~^ and 10~^, respectively. The saturation times in models with rotation are somewhat longer than in models 
without rotation. The clumps formed in the fragmented shell move with velocities that are at least twice as large as 
the escape velocity. The mass of the clumps is ~ 0.1 — 10 Mq, which is lower than the Jeans mass. We conclude that 
the clumps are pressure supported. 

Conclusions. A supernova explosion with energy lO^'^ ergs destructs our model protogalaxy. The clumps formed in the 
fragmented shell are pressure supported. We conclude that protogalaxies with total mass ~ 10^ Mq are unlikely to 
form stars due to high-energy supernova explosions of the first stars. 
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1. Introduction 

The detection of extremely metal-poor stars in our Galaxy 
with the iro n abundance equa l to or less than 10~'^ o f the 
solar value (|Beers et al.L 119921 : IChristlieb erall . l2002f l mo- 
tivated the scientists to put forward pos sible scenarios for 
the fo rmation of such stars. According to lTsuiumoto et all 
(|l999t ) extremely metal-poor (EMP) stars form in a dense 
shell produced by Type II supernova explosions of the first 
stars and accrete metals from the surrounding medium dur- 
ing the subsequent evolution. The formation of EMP stars 
is made possible by fragmentation of the primordial gas in 
a supernova shell due to efficient cooling by molecular hy- 
drogen and HD molecules. Indeed, the formation of molec- 
ular hydroge n increases significantly in the wake of a strong 
shock wave (ISuchkov et al.l . Il983t IShapiro fc Kanej . Il987t) . 
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Moreover, at temperatures lower than 150 K, HD molecules 
become a m uch more efficient coolant than molecular hy- 
drogen (e.g.. IShchekinovl . Il986[ ). 

Subse quent studi es by many authors (see review by 
Nishi fc S usa. 1999; Salvaterra et al.l . l2004HMachida et"aLl . 
20051: [Greif et all . 120071 ) have shown that the likelihood for 
the formation of EMP stars in supernova-driven shells is 
sensitive to both the mass of the dark matter halo (Mh) 
in a primordial galaxy and th e assumed t ype of super- 
nova explosions. For instance, iNishi fc Susa (1999) have 
shown that the shell fragmentation due to Type II super- 
nova explosions (10^^ erg) is possible only in dark mat- 
ter halos with total mass Afh ^ a few x 10^ Mq. Galaxies 
with A/h ^ 10^ Mq may lose its bary onic mat ter before it 
has time to fragment and form stars (jFerraral . [l998). The 
likelihood for this blow-away scenario increases if the first 
stars explode as pair-instability supernovae with the en- 



E. O. Vasiliev, E. I. Vorobyov, and Yu. A. Shchekinov: First supernovae in dwarf protogalaxies 



ergy release of th e order of 10^"^ erg (jBronini et al.l . l2003t 
iGreif et all . l2007l). However, these numerical results were 
confronted bv lMachida et all (J2005l ). who presented a semi- 
analytic model for the evolution of a gas shell produced by 
supernova explosions with energy 10^^ — 10^^ erg. They 
took into account the H2 and HD chemistry and found 
that supernova explosions can induce fragmentation of the 
gas shell and formation of EMP stars in dark matter ha- 
los with Mh = a few x 10^ Mq. The typical mass of the 
fragments in their model is about 1.0 Mf^, w h ich is con- 
sistent with the result by lUehara fc Inutsukal (|2000D and 
similar to the masses o f the observed low-metallicity stars 
(|Christlieb et all . I2002D . 

The likelihood for the formation of EMP stars de- 
pends also on the radial distribution of gas in a primor- 
dial galaxy prior to supernova explosions. This distribu- 
tion can (partly) be det ermined by HII regions aro und the 
first stars. For instance, IWhalen fc Normal] (J2004[ ) consid- 
ered the formation of the HII region in a low-mass dark 
matter halo with mass 10^ Mq and found that the gas is 
efficiently blown away by a supersonic shock wave accoci- 
ated with th e R-type ionizat i on fro nt. Similar results were 
obtained by iKitavama et al.l (|2004l ) for the same mass of 
the dark matter halo. They have also considered more mas- 
sive dark matter halos and concluded that the HII region 
is confined within the virial radius in halos with masses 
^ 10^ MfTi (the ionization front is of the D-type). More re- 
cently, |Kit^^m^^|^^shida (|2005t ) studied numerically the 
destruction of dark matter halos due to the formation of HII 
regions due to radiation of the first stars and subsequent 
supernova explosions. Using a one-dimensional Lagrangean 
hydrodynamics code, they found that the SN shock wave 
remains well inside the virial radius for dark matter ha- 
los with masses larger than 10^ Mq. On the other hand, 
in low-mass dark matter halos with masses of the order of 
10^ Mq the shock wave propagates outside the virial radius. 
These results indicate that the efficiency of halo destruc- 
tion (and, by implication, the efficiency for the formation 
of EMP stars) arc determined not only by the explosion en- 
ergy but also by the gas density distribution and radiative 
feedback from the first stars prior to their explosion. 

It now becomes evident that the formation of EMP stars 
due to fragmentation of supernova driven shells is a com- 
plicated phenomenon that depends on a variety of physi- 
cal conditions in a host protogalaxy, which may vary from 
allowing star formation to shutting it off completely. In 
such circumstances, the construction of increasingly more 
sophisticated numerical models is justified. The majority 
of above mentioned studies are based on semianalytic mod- 
els fe.g.. lSalvaterra et all . I200I : iMachida et al.l . l2005D or on 
one-dimensional numerical h ydrodynamic simulations (e.g., 
IKitavama fc Yoshidal . l2005f ). In the framework of multi- 
dimensional numerical simulations, the priority was often 
given to a careful investigation of dynamical processes lead- 
ing to the destruction of a dwarf protogalaxy by supernova 
explosions, neglecting the chemical processes that control 
the gas thermal properties. In the present paper, we per- 
form axially symmetric numerical hydrodynamics simula- 
tions of high-energy supernova explosions (10^'^ erg) in a 
model dwarf protogalaxy with total (dark matter plus bary- 
onic) mass 1.7 x 10^ Mq. We give a careful treatment to 
the heating and cooling processes and chemical evolution 
of main molecular, ionic, and atomic species. We seek to 
determine the effect of galactic rotation on the dynam- 



ical and chemical evolution of a supernova-driven shell. 
The dark matter halo profile in our models is defined by 
a modifi ed isothe r mal s phere (as observed in local dwarf 
galaxies, iBurkertl. Il995l) . rather than by a cuspy profile 



galaxies, laurkertl. 119 
(|Navarro et all . 119971 ) 



The paper is organized as follows. In Sections 2 and 
3 we describe the model and numerical techniques, re- 
spectively; in Section 4 we discuss several neglected pro- 
cesses; in Section 5 we present results of our numeri- 
cal simulations. The discussion and conclusions are given 
in Sections 6 and 7, respectively. Throughout the pa- 
per we assume a ACDM cosmology with the parameters 
(r2o,^A,^m,^6,/i) = (1-0, 0.76, 0.24, 0.041, 0.73) as in- 
ferred from the Wilkinson Microwave Anisotropy Probe 
(WMAP), and deuterium abundance 2.78 x 10~^, consisten t 
with the most recent measurements (|Spergel et al.l . [20061 ). 



2. Model protogalaxy 

Our model protogalaxy consists of a baryonic component 
surrounded by a spherical dark matter halo. We assume 
that the dark matter halo profile is spherically symmetric 
and is determined by a modified isothermal sphere 
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where rg and po are the core radius and central density, re- 
spectively. Cosmological N-body simulations suggest that 
dark matter halos (upon their formation) have radial pro- 
files that are c uspy in the cen t ral re gions and scale as r^"^ 
at large radii (|Navarro et al.l . Il997l ). On the other hand, 
modeling of the rotation curves in nearby dwarf galaxies 
indicates that the dark matter profil es have a flat c entral 
region and a tail that scales as r~^ (JBurkertl . Il995l ). It is 
poorly known how and when the transition from cuspy 
to flat dark matter halos occurs. Several possible mecha- 
nisms that can facilitate transformation of the profile in- 
clude ran dom velocities of baryons due to a star formation 
feedback (jMashchenko et all . 120061). initial background per- 
turbations of the dark matter (jMikheeva et al.l . l2007t ) and 



ang ular momentum tran sfer from baryo ns to the dar k mat- 
ter (jTonini et al.l . l2006t ). In particular, iMashchenko et al.l 
(|2006[ ) suggest that random bulk motions of gas in small 
primordial galaxies can flatten the central dark matter cusp 
on relatively short timescales, ~ 10* yi'fEJ- Motivated by this 
finding, we assume that dark matter halos in protogalaxies 
at z = 12 have already flattened and attained a configura- 
tion described by equation (P). 

We calculate rg and po as a function of the halo mass 
based on the empirical profile for nearby d warf galaxies 
(IBurkertl . Il995l : ISilich fc Tenorio-Ta glc. 2001') extrapolated 
to the early universe bv lFuiita et al.. (2003. ) 

ro = 8.9 x 10"^ kpc Ml^^h^/^n:;;^^/^{i + z)-\ (2) 

Po = 6.3 x 10^" Mokpc~3 M^^^^h-^^^nmil + zf. (3) 

The dark halo mass (A/h) h^ our numerical simulations is 
set to 10^ Mq, which, at a redshift of z = 12, corresponds to 
3(7 perturbations in the ACDM model for the parameters 
determined from the third year WMAP data. According 



^ For the adopted cosmological model, this timescale corre- 
sponds to 1/4 of the comoving age of the universe at a = 12. 
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to lTegmark et alj ()1997D and lShchekinov fc Vasilie^ (|2006f) . 
protogalaxies with dark matter halos of order 10^ Mq are 
expected to cool down efficiently. The virial radius of our 
model protogalaxy is r^ = 5 20 pc (for relations bet ween 
the virial parameters see e.g., ICiardi fc Ferrarall2004[ ) . 

The initial distribution of the total gas density (pg) is 
found by solving numerically the steady-state momentum 
equations in cylindrical coordinates (z, r) 



<. 


1 dp 


r 


Pgdr 


I dp 
pdz - 


9$h 

dz 



dr 



(4) 
(5) 



where p = pgCg is the gas pressure, Cg = RT/fj, is the sound 
speed, n is the mean molecular weight, ^h is the gravi- 
tational potential of the dark matter halo, and u^ is the 
rotation (azimuthal) velocity of gas. The initial value of fi 
is set to 1.22, which corresponds to a gas of neutral hydro- 
gen and helium with abundances by mass X = 0.76 and 
Y = 0.24, respectively. The gas is assumed to be initially 
isothermal at a virial temperature Tvir = 5900 K. 

In order to solve equations (|4]) and (O for the initial 
gas density distribution, we have to know the initial rota- 
tion velocity of gas, u^q. We determine u^o in two steps. 
First, we calculate the circular velocity of gas Ucirc, which 
is determined exclusively by the gravitational potential of 
the dark matter halo (neglecting the contribution from gas 
pressure gradients) 



dr 



(6) 



Then we introduce a parameter a < 1, which measures the 
relative input of rotation to the total support against grav- 
ity. The rotation velocity is finally determined as m^q — 
aucirc- In the following text, we consider two models: 
model 1 without rotation and model 2 with rotation, for 
which we choose a ~ 0.4. The resulted total rotational 
energy inside the virial radius is Sj-ot = / w?g p dV w 
0.31 X 10^^ erg. The corresponding ratio /3 of rotational to 
gravitational energy is equal to 0.17. We note that /3 ~ a^ 
because -Erot/^-gr — '"^/''^circi where Egi- is the gravita- 
tional energy. We note that our model galaxy has a modest 
amount of rotation and most support against gravity comes 
from gas pressure gradients. 

The numerical procedure for solving the steady-state 
equations ^ and ([5]) to obtain the equil i brium gas density 
distribution is given in IVorobvov et al.l (|2004f ). There is, 
however, one important difference. The total masses of gas 
and dark matter halo in protogalaxies are linked through 
the relation Afg = AI\^{ni,/flm)- In order to satisfy this 
relation, we iterate solutions until the resulted gas density 
distribution has attained the mass of about 1.7 x 10^ Mq 
inside the virial radius. 

Once the equilibrium gas density distribution is con- 
structed, we let our model galaxies evolve out of equilib- 
rium. The inner densest regions are characterized by largest 
cooling rates, and they cool and contract on much shorter 
time scales than the outer regions. We stop this process 
when the temperature in the center drops below 500 K. cor- 
respondingly the density increases about 50 times in com- 
parison with the central equilibrium value. This happens at 
t > 1.0 Myr in model 1 and t > 1.4 Myr in model 2. The 



resulted distributions of gas density (top panel), tempera- 
ture (middle panel), and infall velocity (bottom panel) in 
model 1 (solid line) and model 2 (dashed and dotted lines) 
are shown in Figure [T] More specifically, the dashed and 
dotted lines show the radial and vertical profiles in model 2, 
respectively. We note that in the non-rotating model 1 the 
distributions are identical in all directions. 

The infall velocities in Figure [T] indicate that only the 
inner regions of our model galaxy r ^ 50 pc are driven 
out of equilibrium, while the rest of the galaxy is intact. 
This out-of-equilibrium region remains confined to the in- 
ner 100 pc, even if we extend our numerical integration for 
another several million years. As a result, the radial gas 
density profiles in model 1 and 2 are distinct throughout 
the bulk of the galaxy. In particular, the radial gas density 
profile in model 2 (with rotation) is shallower (oc r~^'^) 
than in model 1 without rotation (ex r~^'^), except for the 
innermost region. The difference in the central gas density 
between the models with and without rotation can reach a 
factor of three. However, the masses enclosed inside the 
virial radius (w 520 pc) are approximately the same in 
both models. The vertical gas distribution in model 2 (dot- 
ted line) is characterized by progressively smaller values 
(at same radii) than in model 1 (solid line). This implies 
that supernovae remnants in rotating protogalaxies are ex- 
pected to expand to a larger distance in the vertical direc- 
tion, which may result in a partial loss of gas by a parent 
protogalaxy. 

Once the central gas temperature drops below 500 K, 
we release 10^^ erg of thermal energy in the central sphere 
with radius 5 pc. Such energetic supernovae are expected to 
result fror n the pair-instabili t y exp losion of massive metal- 
free stars (iHeger fc Wooslevl . |2002| ). This approach is sup- 
posed to reasonably mimic the formation of a first star and 
its subsequent explosion. The gas temperature in the cen- 
tral sphere is set to 2.5 x 10"*^ K, typical for HII regions, and 
gas density is adjusted to maintain the pressure balance be- 
tween the interior and immediate exterior of the sphere. We 
assume that this sphere was formed by a collective action 
of stellar wind and ultraviolet radiation from a massive star 
(the supernova progenitor) that had resided in the galactic 
center before the explosion (we skip this phase in our nu- 
merical simulations) . Since the physical size of a shock front 
at initial stages is considerably smaller than our numerical 
resolution (~ 1 pc). we allow for the gas to initially evolve 
adiabatically and switch on cooling only after 6 x 10^ yr 
of evolution. In what follow we will loosely define this hot 
sphere as a supernova ejecta, although it differs physically 
from the ejecta in realistic supernovae explosions. The ini- 
tial composition of species inside the sphere corresponds to 
a fully ionized gas. 



3. Numerical model 

3.1. Main equations 

The dynamics of the gaseous component in our model pro- 
togalaxy is followed by numerically solving a usual set of 
hydrodynamic equations in cylindrical coordinates (z,r, 0) 



9pg 
dt 



V • (pgw) = 0, 



di 



(pgu) + V • {pgu ■ u) = -Vp - PgV$h, 



(7) 
(8) 
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Fig. 1. Profiles for the gas density (top panel), tempera- 
ture (middle panel), and infall velocity (bottom panel) in 
the non-rotation model 1 (solid line) and rotating model 2 
(dashed and dotted lines). In particular, the dashed and 
dotted lines show the radial and vertical profiles in model 2. 
respectively. 



de 

dt 



— + V-(eM) = -p(V 



A, 



(9) 



where u — zu^ + rur + 4>u^ is the gas velocity in cylindrical 
coordinates, V = zd/dz -I- rd/dr is the spatial derivative 
in cylindrical coordinates modified for the adopted axial 
symmetry, A is the cooling rate, and e is the internal energy 
density. The hydrodynamic equations are closed by an ideal 
equation of state p = (7 — l)e, where 7 is set to 5/3. A 



(small) amount of numerical viscosity is added to the code 
(according to the usual prescription of von Neumann and 
Richtmycr) in order to smooth out shock and discontinuities 
over two computational zones. 

Equations ([7][9]) arc solved using a finite-difference, 
operator-split code, which ap plies the cons i stent trans- 
port technique introduced by iNorman et al.l (|1980[ ). The 
code performs well on the standard test suite including 
the Sedov point explosion test. The advection is treated 
using a third-order picccwise parabolic advection scheme 
(Collcla & Woodward, 198j). This scheme is known to be 
more accurate in handling the shocks and contact discon- 
tinuities than a commonly used v an Leer advecion scheme 
(see e.g.. IStone fc Normanl . 11993 ). The computational do- 
main has a size of 750 pc in both the vertical (z) and hor- 
izontal (r) directions. Cosmological effects (expansion and 
bias) are not expected to significantly influence the results 
on such spatial scales. We assume an equatorial symme- 
try for simplicity. The numerical resolution is 780 x 780 
grid zones, which corresponds to a spatial resolution of 
« 1.0 pc in both coordinate directions. In our simulations 
self-gravity of the gas is neglected and the dark matter po- 
tential is fixed. In a few test runs, we included the self- 
gravity of the gas but found it dynamically unimportant. 

3.2. Cooling processes and chemical reactions 

The gas component of our model protogalaxy consists of a 
standard set of species: H, He, 11+ , H", II2, 112", D, D+, 
and HD. We also make use of the charge conservation law 
for electrons. The initial number densities of neutral hydro- 
gen and helium relative to the total number density of gas 
rig = pg/{fxmn) are a;[H] = 0.9999 x 0.76 and a;[He] = 0.24, 
respectively. The initial number densities of ionized hydro- 
gen, molecular hydrogen, and HD molecules relative to n^ 
are a;[H+] = lO"", xpz] = 10"^ and x[HD] = lO"!", re- 
spectively. The initial relative number density of deuterium 
is equal to the cosmological value. The relative number den- 
sities of the other species in our sample are set initially to 
a negligible value. 

The cooling rates are computed separately for temper- 
atures below and above 2 x lO''' K. In the low-temperature 
regime, the cooling rate A in the energy equation ^ in- 
cludes typical radiative losses in the primordial plasma. 
Cooling due to recombination an d col li siona. 1 excitation of 
atomic hyd rogen is taken from ICenI (|1992f ) . It is worth 
mentioning iGalli fc Pallal ()1998f ) H2 cooling function al- 
most concides with the total cooling functi on including H- 
H2 a nd e~-H2 collisions for Xe. < 10~^ () Glover fc Abell . 
|1998[ ). which is the case in our sim ulation. Cooli n g due 
to molecular hydrogen is taken from iGalli fc Pallal ( 19981 ) 
and modified for the t emperature regime near the CMB 
radiation t emperature ([Varshalovich fc Khersonskiil 1X9761 : 
iPuv et all . Il993t iLe Bourlot et all . Il999t) . Coohng due to 
H D molec u les is comp uted according t he pre scription given 
in iFlowed (|200Gf ) and iLipovka et all ^B)- In the high- 
temperature re gime, the cooling rates for zero metallicity 
are taken from [Sutherland fc Dopital (|1993( ). 

We assume that our atomic, ionic, and molecular species 
are collisionally coupled to each other and share a common 
velocity field, which eliminates the need for solving separate 
equations of motion for each species. Hence, in addition to 
the usual hydrodynamic equations ([7]Oi we have to solve 
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Table 1. Chemical reaction rates 



reaction 




Reference 


H+ + e~ - 


->H + h!y 


Galh & Palla (1998) 


H + e- -^ 


H-+h// 


- 


H- + hv - 


-^H + e" 


- 


H^ +H- 


+ Ha +e" 


- 


H-+H+ 


^ H+ + e" 


- 


H-+H+ 


^2H 


- 


H+ + H- 


+ H+ + hi/ 


- 


H+ + H- 


+ Ha + H+ 


- 


H++e- - 


-^2H 


- 


Ha + H+ - 


-^H++H 


- 


Ha+c- - 


^ 2H + e" 


- 


Ha + H^ 


■3H 


MacLow & ShuU fl986) 


H + c" ^ 


H+ + 2e- 


Abel et al. fl997l 


H^ +0" - 


-^ H + 2e- 


- 


H^ +H- 


* 2H + e" 


ShaDiro & Kane f 19871 


D+ + H2 


~> HD + H+ 


Galh & Palla (1998) 


H++HD 


-> Ha + D+ 


- 



the continuity and rate equations for the mass densities (pi) 
of each of the species 



dt 



V • {piU) ^ fcjksPjPk - /^kiPkPi, 



(10) 



where the right-hand terms are the sources and sinks due 
to chemical reactions. The list of reactions and references 
to the corresponding formation (fcji) and destruction (fcki) 
rates are given in Table [1] We note that wc do not solve 
equations (fTO|) for D+ and D, since the ratio of ionized 
versus neutral deuterium densities can be implicitly derived 
from the ratio of ionized vers us atomic hydrogen densities 
(jVasiliev fc Shchekinovl . l2003l ) . We have also neglected the 
helium kinetics. 

The hydrodynamic equations ^]^ and the chemical re- 
action network (jlOp are numerically coupled using the fol- 
lowing strategy. First, the cooling rate A is computed and 
the global hydrodynamic time step Aih is derived as 



1 



1 



At 



CFL 



At 



cool 



(11) 



where Ai^pL is a usual hydrodynamic time step due to the 
Courant-Friedrichs-Lewy condition and Aicooi = O.le/A is 
a characteristic cooling time. Note that Ath should be min- 
imized over all computational cells. 

Once the global hydrodynamic time step is calculated, 
the solution of equations ([7|)- pII)l proceeds as follows. First, 
we solve the continuity and momentum equations for the 
total gas density pg (cqs. ^ and ([5]), respectively) using 
a technique described in Section 13.11 The solution of the 
chemical reaction network (jlO|) and energy balance equa- 
tion ^ is split into two parts. The update of internal energy 
density (e) and mass densities of each species (pi) due to 
advection is done by solving the following continuity equa- 
tions 



^ + V • {p,u) = 0, 



(12) 
(13) 



The solution technique is exactly the same as for equa- 
tion ©. 



The remaining update of pi due to chemical reactions 
is combined with the update of e due to cooling and com- 
pressional heating 



dpi 
dt 
de 
di 



fcjkPjPk - fckiPkPi 

-p(V • It) - A. 



(14) 
(15) 



To improve the accuracy, we use a predictor-corrector 
scheme. At each step of the predictor-corrector method, 
equations ([T^ are integrated using a fifth-order Runge- 
Kutta-Carsp meth od with the adaptive step size control 
([Press et al.l . 11993 ). The global integration time step is de- 
termined by Aih- We note that in the high-temperature 
regime (T > 2 x 10^ K) the mass densities of atomic and 
ionized hydrogen arc computed adopting a local thermal 
equilibrium. More specifically, the ionization and recombi- 
nation rates for hydrogen are set equal. The mass densities 
of the other four species (H~, H2, H^, and HD) arc set to 
a negligible value. 



4. Neglected processes 

In our numerical simulations, we have neglected several 
physical processes that may be important in some dwarf 
protogalaxies. For instance, a massive star - a likely su- 
pernova progenitor - emits enormous ultraviolet (UV) flux 
that ionizes the surrounding gas and forms the HII region. 
The radiative pressure of UV photons from massive stars 
is so strong that it can blow the gas away from the cen- 
ter of a low-mass (^ 10^ Mq at z = 20) protogalaxy 
(IWhalen fc Normanl . |2004 ). Our model protogalaxy has a 
more massive dark matter halo, M = 10^ Mq, but is lo- 
cated at a lower redshift, z = 12. An HII region formed 
inside such a protogalaxy is expected to be confined within 
10 pc. This estimate is obtained by assuming the emission 
rate of ionizing photons from the first stars to be lO^'' s~-^ 
(Schaercr, 200^) and background density n = 10 cm""^. It 
is worth noting however that the size of HII region depends 
on the dark matter profile of a protogalaxy: it is expected 
to be larger in protogalaxies with a flat density profiles, and 
lower in the cuspy profiles. Nevertheless, in order to avoid 
an overestimate of the effects of the HII region, we set its 
size in our numerical simulations to 5.0 pc. 

Shock waves, a nother source of ionizing radiation 
(IShuU fc Silkl . Il979f ). are efficient only in the adiabatic 
phase of expansion, when the temperature is much higher 
than 10^ K. They cannot infiuence significantly the forma- 
tion of molecules at later times. 

Another feedback process that can affect the chem- 
istry of gas is the radiation in the (11.18 — 13.6) eV band, 
which photo-dissociates H2 molecules. Although a photo- 
dissociative region creat ed by the rad i ation of a massive 
star can be quite large (jHaiman et al.l . Il997f ). the typical 
relative molecular hydrogen density in this region is quite 
small, about 10~^ (taken as an initial value in our numerical 
simulations). Obviously, the primary source of dissociative 
photons ceases to exist after the supernova explosion and 
we expect that the production of dissociative photons by 
the hot gas of a bubble created by the explosion is negli- 
gible. Moreover, H2 molecules cannot be photo-dissociated 
inside cold and dense regions o f a supernova-driven she ll 
due to the self-shielding effect (jDraine fc Bertoldil . |1996[ ). 
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All this argues against photo-dissociation as an important 
mechanism for H2 destruction after the supernova explo- 
sion. 

A positive feedback on the abundance of H2 molecules 
can come f rom the X-ray photons produ c ed by supernova 
explosions (JFerraral Il998t iHaiman et all . I1997D . However, 
other coacting mechanisms are likely to be more important. 
For instance, it is well established that the H2 formation is 
very efficient behind strong shock waves (jSuchkov et al.l . 
ll983HShapiro fc Kan'i . ll987n . As a result, the relative den- 
sity of H2 molecules increases very rapidly below 10'' K 
and saturates at a so-called "universal" value of about 
a:[H2l ~ (1 - 2) x lO''^ at temperature ~ 3 x 10^ K 
(jOh fc Haimanl . l2002f) . Thus, the X-ray photons have only a 
minor influence on the H2 abundance in a supernova-driven 
shell but they can strongly influence the H2 abundance in 
the gas ahead the shell. However, this is not expected to 
change the picture as a whole because shock waves created 
by superovae with the energy release of order 10^^ erg are 
strong. 

To summarize, we do not expect that these processes are 
significant for our model protogalaxy. In a future work, we 
plan to include radiative transfer in the numerical hydrody- 
namics simulations to check the validity of our assumptions. 



5. Results 

In this section we perform a comparative study of the long- 
term evolution of supernova remnants in the non-rotating 
and mildly rotating galaxies. We note that the equilibrium 
configuration of gas is not preserved exactly during our nu- 
merical simulations. A gradual drift from the initial equi- 
librium (though very small) introduces a seed perturbation 
to the gas density, which is necessary for the instabilities 
to grow. These perturbations can be regarded as initial gas 
inhomogeneities that are always present in protogalaxies 
before supernovae explosions. 

Figure 2 presents snapshots of the gas density distri- 
bution in model 1 (top row) and model 2 (bottom row) at 
four consecutive times after the supernova explosion. In the 
non-rotating model 1, a thin shell of compressed gas sepa- 
rating the hot supernova ejecta from the unperturbed gas is 
clearly seen at i = 1 Myr. The shock wave has decelerated 
by this time and its position coincides with the position 
of the shell. On the contrary, in the rotating model 2 the 
shock wave lies ahead of the shell at the same evolution- 
ary time. This is a consequence of a lower gas density in 
the center of model 2 (see dashed and dot-dashed lines in 
Fig- HI) ■ It takes a longer time for the shock wave to decel- 
erate and merge with the shell in model 2. The shock front 
in model 2 is profoundly elongated along the vertical axis 
due to a non-spherical initial distribution of gas. 

When the characteristic cooling time becomes shorter 
than the dynamical time (the age of a supernova remnant), 
an expanding shell becomes unstable to th e Rayle igh- Taylor 
and Kelvin-Helmholtz instabilities (jGulll . Il973l ). As a re- 
sult, small ripples that distort a spherical shape of the shell 
appear in model 1 at f = 1 Myr. Moreover, because the 
shell is radiative, we expect that the thermal instability 
plays a non- negligible role in the formation of the ripples. 
The subsequent evolution of the shell is governed by the 
Rayleigh- Taylor instability, which acts mostly in the com- 
pressed gas of the shell outside the interface between hot 
supernova ejecta and the shell of compressed material. The 



ejecta itself experiences only large scale distortions, which 
seem to be inefficient in mixing the metals throughout the 
shell. The characteristic time for the development of the 
Rayleigh- Taylor instability is shorter for steeper initial gas 
density profiles and vice versa. Figure [1] indicates that both 
the gas density distribution in model 1 and the vertical gas 
density distribution in model 2 have profiles similar to r~^, 
where r* = (r^ -|- z^)^/^ is the distance from the galactic 
center. On the other hand, the radial gas density distribu- 
tion in model 2 is noticeably shallower and follows an r~^-^ 
profile. Hence, we expect the Rayleigh- Taylor instability 
to grow faster in the non-rotating model 1. This is indeed 
seen in the top row of Fig. [5] - the shell has lost its spher- 
ical shape by t = 4 Myr and prominent spurs (or fingers) 
start to grow into the unperturbed medium. Model 2 shows 
little spurs at the same evolutionary time, though the shell 
has already started to show first signs of instability. In the 
end of numerical simulations, model 1 develops consider- 
ably longer spurs than model 2. However, the number of 
the spurs is smaller in model 1 due to a larger characteris- 
tic length of instability. We note that the Rayleigh- Taylor 
instability is expected to grow faster in protogalaxies with 
cuspy dark matter halos due to a steeper gas density pro- 
file than in protogalaxies with halos characterized by a fiat 
central region. However, this effect may be p r esent only 
in massive protogalaxies (JKitavama fc Yoshidal . |2005[ ). be- 
cause steep gas density profiles in low-mass galaxies with 
cuspy dark matter halos are likely to be destroyed by ion- 
izing radiation from the first stars. 

Figure [2] shows that a,t t = 12 Myr some spurs are 
found outside the virial radius (520 pc) in the both models. 
This implies that a fraction of the baryonic mass is blown 
away by the supernova explosion. To calculate this frac- 
tion (/out), we notice that the spurs are characterized by 
a systematically lower temperature than the unperturbed 
medium. We use this property and find /out by summing 
up the gas mass in the computational cells occupied by the 
spurs at t = 12 Myr. We count only those computational 
cells that lie outside the virial radius. As the spurs cross the 
virial radius, they sweep up some of the pristine gas. Since 
we are only interested in the blown-away gas, we subtract 
the input from this unperturbed gas. The resulted fraction 
/out as a function of time is shown in Fig. [8] by the solid 
(model 1) and dashed (model 2) lines. It is evident that the 
rotating protogalaxy (model 2) losses roughly twice as much 
baryonic mass as the non-rotating protogalaxy (model 1). 
This is a consequence of a shallower initial gas density pro- 
file in the rotating protogalaxy. 

The spurs have a complicated internal structure. The 
central regions (or cores) of the spurs are characterized 
by temperatures lower than 1.5 x 10"^ K and densities 
roughly ten times larger than those of the neighbouring 
unperturbed gas. Because the spurs move almost ballisti- 
cally through the unperturbed gas, bow shocks form around 
them, creating envelopes of shocked gas. The envelopes are 
hotter than the cores and are characterized by tempera- 
tures about 10^ K and densities roughly four times larger 
than those of the neighbouring unperturbed gas. 

Is the blown-away gas lost to the parent galaxy? 
Figure [4] shows the gas velocity field superimposed on the 
gas density distribution at t = 12 Myr. We find that the 
spurs are characterized by mean mass-weighted velocities 
of the order of 26 km s~^ in model 1 and 22 km s^^ in 
model 2. In the latter model, higher velocities are found 
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Fig. 2. Logarithmic distribution of density (in cm^'^) at i = 1, 4, 8, 12 Myr after SN explosion with energy Esn 
Xo53 gj-g JQ j^aio with mass M = 10^ Mq for the spin parameter (3 = 0- upper and /3 = 0.17 - low row of panels. 



near the vertical axis but it may be a numerical artifact 
due to the axisymmetric nature of our numerical simula- 
tions. The spurs appear to have already escaped the galaxy. 
Indeed, the escape velocity at the virial radius in both mod- 
els is We = (25r'^v)^^^ ~ 13 km s~^, where g^ = GMy{ry)/rl 
is the radial gravitational acceleration at the virial radius 
and My{r^) is the halo mass inside r^. It is obvious that the 
mean mass-weighted velocity of the spurs is at least twice 
as large as the escape velocity at the virial radius, which 
implies that the spurs might have escaped our protogalaxy. 

Figure shows the distribution of gas temperature 
(left column), relative number densities of molecular hy- 
drogen (middle column) and HD molecules (right column) 
in the end of numerical simulations at i = 12 Myr. The 
top/bottom rows correspond to model 1 /model 2. It is ev- 
ident that low gas temperatures (below 10^ K) are found 
in the shell and the spurs, where cooling takes place due 
to II2 and HD molecules. In particular, the lowest tem- 
peratures found in the spur cores are of order 500 K. The 
relative number densities of molecular hydrogen and HD 
molecules in the spur cores are quite large, approximately 
10~^ and 10~^, respectively. We note that the quoted H2 
relativ e number densi t y is c lose to a so-called "freeze-out" 
value (|0h fc Haimanl . |2002[ ). A sharp decrease in the H2 
relative number density in the envelopes is explained by 
efficient destruction of H2 molecules due to collisions with 
hydrogen atoms via the following reaction: H2 -I- H — > 3H. 
Shock waves produced by supersonic spurs heat the com- 
pressed gas to temperatures above 10^ K, at which the col- 
lisional destruction of molecular hydrogen becomes domi- 
nant. There is no significant difference in the temperatures 
and relative number densities of molecular species found 
between the two models. That is not unexpected because 
strong shock waves lead to similar final temperatures and 
relative number densities of species in the gas. 




t, Myr 



Fig. 3. Fraction of the baryonic mass blown away by a su- 
pernova energy release of 10^^ ergs in model 1 (solid) and 
in model 2 (dash). 



Figures [6] and [7] present the temporal evolution of differ- 
ent molecular hydrogen tracers in model 1 and model 2, re- 
spectively. In particular, filled squares show the total molec- 
ular hydrogen mass Mh2, the filled and open circles yield 
the total mass of gas with the relative number density of 
molecular hydrogen a;[H2] > 5 x 10""' and x[H2] > 10"'^, 
respectively. In addition, we keep track of the gas mass with 
temperature < 10'^ K (filled triangles) and < 500 K (open 
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Fig. 4. Velocity field overlaid on the logarithmic density map. Left panel: protogalaxy without rotation, right: 
rotation. The arrow in the upper-right corner corresponds to the value of velocity 50 km s^^. 
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Fig. 5. Logarithmic distribution of temperature {left panel), H2 (middle) and HD (right) abundances at time t = 12 Myr 
after the explosion of SN with energy Esn = 10^^ erg in halo with mass M = 10^ Mq and the parameter /3 = (upper 
panels) and (3 = 0.4 (lower panels). 



trinagles). All molecular hydrogen tracers are calculated in- 
side the computational domain. 

The comparison of Figs. [H] and [7] shows that the H2 
traces saturate during the evolution. The saturation is ex- 
plained by the fact that we consider the gas evolution be- 
hind strong shock waves, where the H2 relative density 



rapidly reached a "universal" value ^ 10"'^ (see section 4). 
The saturation times in model 2 are systematically longer 
than in model 1. This can be attributed to longer cool- 
ing times in model 2 due to a shallower initial gas den- 
sity profile. The filled circles in Figs. [6] and [7] indicate that 
the total mass of gas with x[ii2] > 5 x 10""' saturates at 
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Fig. 6. The total mass of gas in computational domain for 
model 1 with H2 abundance higher than x[H2] = 5 x lO^'' 
(filled circles), a;[H2] = 10"'^ (open circles), with tempera- 
ture lower than T < 10^ K (filled triangles), T < 500 K 
(open triagles). The line with filled squares represents the 
total mass of molecular hydrogen. The dashed line shows 
the gas mass contained in fragments with density log n > 
-0.25 and temperature T < 5 x 10^ K. 



(2 — 3) X 10^ Mq, which corresponds to ~ 10% of the total 
baryonic mass in our model galaxy. The total gas mass with 
a^[H2] > 10~^ (open circles) saturates at ~ 10^ Mq. The gas 
mass with temperatures below 500 K approaches an upper 
limit of ~ 2 X 10^ Mq. The total mass of molecular gas sat- 
urates at approxim ately 300 — 400 A/0, which is close the 
estimates made by iFerraral (|1998D . It is also evident that 



the H2 tracers in model 2 are characterized by saturated 
values that are systematically larger than in model 1. This 
is due to the fact that the gas mass that crosses the shock 
wave front is larger in model 2 than in model 1. 



6. Discussion 

We have considered a SN explosion in protogalaxies with 
and without rotation, and found that rotation changes ge- 
ometry and dynamics of the flow after the explosion, but 
does not affect significantly gas chemistry of the final state. 
Some differences can be found in temporal evolution, how- 
ever the saturated values of molecular abundances in all 
considered models in the two cases lie in a very narrow 
range. It is quite expected that protogalaxies with the ra- 
tio of rotational to gravitational energy larger than that 
in model 2 (/3 = 0.17) are strongly violated in the ver- 
tical direction by radiation from the progenitor of a SN. 
This is basically connected with the fact that the cen- 
tral gas density decreases with the parameter /3 approxi- 
mately as ric ^ exp(— /3), and the characteristic size of a 

photo-ionizationally violated region varies as R ^ Uc . 



Fig. 7. The same as in Fig. 6 for model 2. 

Moreover, disks in such halos are more extended in the ra- 
dial direction with a flat radial density proflle 



n ^ r , 

where A ~ —2/3 + 2, versus 
-,-2 



(16) 



(17) 



in the vertical direction. This means that photoionizing ra- 
diation propagates and evacuates gas mostly perpendicular 
to the plane around the galactic center. As a result, when 
a SN explodes, its shell will expand through this cave, such 
that a significant fraction of gas can be blown away into 
the intergalactic medium. One can therefore speculatively 
conclude that the effects from ionizing photons of the pro- 
genitor of a SN will be more pronounced in rotating galaxies 
in the form of jet-like outfiows. 

6.1. Evolution of fragments 

One can see from Figs. [5] and [5] that the typical radial 
length of the most cold and dense regions is several par- 
sees. The density of such clumps is about an order of mag- 
nitude higher than that of the surrounding gas, whereas 
the temperature is about an order lower, thus the clumps 
are in pressure equilibrium with the environment. The 
clumps form due to dcsintegration of the shock wave under 
Rayleigh- Taylor instability which develops when gas behind 
the shock front starts cooling rapidly and the front decel- 
erates. A typical size of fragments is expected to be close 
to the thickness of the compresses gas behind the front at 
the moment, when it becomes unstable. 

Among possible mechanisms of cloud destruction strip- 
ping of the external layers of clouds seems the most effi- 
cient under the conditions of interest. This pro cess operates 
mostly by Kclvin-Helmholtz instability (e.g., iKlein et al.l 
,1994|) continuing on shorter scales of a single cloud. The 
typical stripping time is therefore of the order of Kelvin- 



Helmholtz instability t 
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Il961t) . where x is the density ratio between the cloud 
and the background medium. For typical parameters v = 
20 km s~^, X ~ 10 and the wavelength of the order of the 
clump radius ao ~ 10 pc, the destruction time is ^ 3 Myr. 
This is short compared to the dynamical time, and from this 
point of view dense clumps should be destroyed quickly. 
However, the radiative cooling time is of the same order 
tc ~ 1 — 3 Myr, which means that the density increase al- 
ways connected with the radiative cooling can inhibite the 
destruction through stripping, so that the clumps can sur- 
vive on longer dynamical time. The clumps however asymp- 
totically are destroyed, which is seen from the fact that 
the mass contained in relatively dense (n > 0.56 cm""^) 
and cold (T < 5 x 10'^ K) fragments decreases at f > 3 
Myr as shown in Fig. 6 and 7 by dashed lines. Moreover, 
a typical mass of the most dense (n > 1.8 cm"-^) clumps 
M -- 0.1 - 10 Mq (see Fig. 8) is much smaller than the 
Jeans mass for the corresponding conditions (n ^ 1.8 cm""^ 
and T ~ 0.5 — 1 X 10'^ K). All this means that protostellar 
clouds do not form in the shell unless the clumps merge. 

In this connection it is worth mentioning the conclusion 
bv lSalvaterra et all (|2004D that the shell formed by a SN in 
the conditions similar to that considered here can fragment 
only through gravitational instability on to larger masses. 
One can think that the two instabilities - Kelvin-Helmholtz 
and gravitational - complement each other, in the sense 
that while Kelvin-Helmholtz instability of an expanding SN 
shell breakes it on to separate fragments (clouds), the grav- 
itational instability stimulates merging of these low-mass 
fragments, and essentially can result in formation of pro- 
tostellar clouds of sufficiently large masses. This possibility 
depends, however, on interelation between the characteris- 
tic times of gravitational instability and destruction of the 
fragments through stripping. As seen from Fig. 4 the rela- 
tive velocities between the fragments and the ambient gas 
remain asymptotically quite high, 20-30 km s~^, which re- 
sults in a relatively short stripping destruction time: ~ 1 
Myr ~ much shorter than Jeans time. From this point of 
view the process of star formation in an expanding SN shell 
meets difficulties: from one side, low-mass fragments origi- 
nated through disruption of the shell are sub- Jeans and can- 
not form protostellar condensation, on the other side, de- 
struction of these fragments by stripping prevents them to 
be gravitationally collected into massive overcritical clouds. 
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Fig. 8. Number of fragments in the expanding shell at 
i = 12 Myr for a protogalaxy without (left) and with (right) 
rotation; panels from the uppermost to the lowermost cor- 
respond to the fragments with log n > —0.25, 0, 0.25, 
respectively; the mass of the fragments was calculated as 
pS^''^, where S is the area of a filamentary fragment on 
the "radius-height" plane, p is the mean density in it; note, 
that the fragments of low density do look more filamentary 
and irregularly shaped than the dense ones, what therefore 
explains a higher spread of masses of fragments with lower 
density limit. 



6.2. Mixing of metals 

Supernovae are recognized as a factory of metals and 
dust, although many issues related to the efficiency of 
metal mixing and production of dust particles inside the 
ejecta, destruct ion of dust in supernovae shells are still un- 
der discussion (iMadau et all. 12001 1 Suger man et al.l . l2006t 
iMeikle et al.l . l2007l : IVenkatesan et al.. . .2006l) . In particular, 
two conflicting conclusions about t he efficiency of dust pro - 
duction in SN 2003g d event by ( ISugerm an et all |2006[ ) 
from one side and by (jMeikle et al .. 2007) from the other 
have to be mentioned. As far as mixing of metals in a SN 
driven shell is concerned, we believe that metals mostly 
confined in the border between the hot bubble and the 
surrounding interstellar gas, do penetrate and mix in the 
swept-up shell very slowly. Mixing of metals in these con- 
ditions is deter mined by Rayleigh - Taylo r instability in the 
interface layer (jWang fc Chevalien . l200lh which develops on 



times longer than the age of the remnant. Indeed, for de- 
celerating shells the criterion for Rayleigh- Taylor instability 
reads k < 1/6 A R, where Ai^ is the thickness of t he sh ell, k 
iswavenumber (jVishniad . ll983trVishniac fc Rvul . ll989t ). On 
the other hand, the instability increment is 7r,t ~ A/A:|a|, 
where a is the deceleration. This gives for a radiative rem- 
nant of the age t the restriction 7RTi < v2, which means 
that the perturbations amplitude can grow at most by fac- 
tor of e^^ (JKas'vanova fc Shchekinovl ■ l2005^ ■ and therefore, 
mixing of metals through the whole volume of the swept- 
up gas seems to be rather inefficient. From the point of 
view of our simulations this means that the effects of met- 
als on cooling and chemistry of the swept-up gas can be 
neglected in the initial expansion t ^ 1 — 3 Myr. Later on, 
at f >, 3 Myr, the shell enters an accelerating phase so that 
Rayleigh- Taylor instability grows faster. However, the in- 
terface between the shell and the metal enriched ejecta re- 
mains weakly distorted on large spatial scales, and therefore 
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the metals remain mostly locked in a relatively narrow layer 
around the interface. 



7. Conclusions 

In this paper we have considered numerically the effect of 
energetic supernovae explosions (10^'^ erg) in non-rotating 
and rotating protogalaxies with the total mass lO'' Mq at 
a redshift of z = 12. The assumption of axial symmetry 
allowed us to evolve the supernova driven shell for several 
tens of million years. Specifically, we find the following. 

— The process of the shell destruction is different for non- 
rotating and rotating protogalaxies. Fingers or spurs 
formed due to various instabilities (mainly, due to the 
Rayleigh- Taylor instability) grow faster and can be 
found on larger distances in non-rotating galaxies. That 
is explained by a steeper initial gas density profile in 
non-rotating protogalaxies and, as a consequence, by a 
faster characteristic time-scale for the development of 
the Rayleigh- Taylor instability. 

— The supernova evacuates a significant portion of the ini- 
tial baryonic mass to radial distances comparable with 
or larger than the virial radius. At i = 12 Myr after 
the supernova explosion some of the fragments of the 
destructed shell have attained mean mass-weighted ve- 
locities of the order of 30 km s^^ and 25 km s~^ in non- 
rotating and rotating protogalaxies, respectively. These 
values are at least twice as large as the escape velocity 
at the virial radius, which implies that the fragments 
eventually may escape the protogalaxy. We estimate the 
fraction of this blown-away gas in the end of our nu- 
merical simulations to be about 10% and 20% of the 
initial gas mass for models without and with rotation, 
respectively. The re-collapse phase of the shell sets in 
faster in models without rotation than in models with 
rotation. These results are mostly consistent with the 
intuitive expectations as far as the difference between 
models with and without rotation is concerned. 

— The relative number densities of molecular hydrogen 
and HD molecules in the fingers and spurs are found 
to be quite large, approximately 10~^ and 10~^, respec- 
tively. The typical temperature in the spur cores is of 
the order of 500 K at t <L 8 Myr after the supernova 
explosion. 

— The total gas mass with a;[H2] > 10~^ saturates at 
2 X 10^ 7\/0, which corresponds to approximately 10% 
of the total baryonic mass in our model galaxy. The cor- 
responding total mass of molecular hydrogen is about 
300 — 400 Mq . The saturation time-scales in the model 
with rotation are systematically longer than in the 
model without rotation. This can be attributed to longer 
characteristic cooling times in the model with rotation 
due to its shallower initial gas density profile. 

Finally, we would like to note that the typical masses 
of most fragments (assuming their spherical symmetry) are 
~ 0.1—10 Mq. For the typical densities and temperatures in 
the fragments to be ~ 0.5 — 1 cm~^ and (0.5 — 1) x 10'^ K, 
respectively, these masses are strongly sub- Jeans and the 
fragments are expected to be pressure-supported. Their fur- 
ther evolution depends on both the efficiency of cooling 
and destruction due to the Kelvin-Helmholtz instability. 
We do not expect that the low-mass stars can be formed 



in such conditions as suggested by MacKey et al. 2003, 
and Salvaterra et al. 2004. In our opinion, a more feasi- 
ble mechanism for low-mass, metal-poor star formation is 
related with the re-collapse of a supernova bubble in proto- 
galaxies, whose total gravitational binding energy is much 
larger than supernova energy. 
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